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Abstract 

We calculate time correlation functions in the Bak-Sneppen model 
(Phys. Rev. Lett. 71 4083 (1993)), a model showing self-organised 
criticality. For a random neighbour version of the model, analytical 
results are presented, while on a one dimensional lattice we give nu- 
merical results. The power spectrum of these correlation functions 
shows 1//- behaviour in both cases. 
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A large diversity of physical systems show 1 / /-noise . The power spectra 
of time correlation functions of such systems show powerlaw behaviour /"^ 
over many orders of magnitude with an exponent (3 in the range .6 ~ 1.6. 
A possible explanation for the wide occurence of this phenomenon was put 
forward in a paper entitled "Self-organised criticality: An explanation of 
1// noise" p|. In that paper, Bak, Tang and Wiesenfeld argue that 
many open non-linear dynamical systems with large number of degrees of 
freedom evolve to a state where they show critical behaviour characterised 
by powerlaw correlations both in space and time. Bak,Tang and Wiesenfeld 
(BTW) illustrate their ideas using a simple model, the so called sandpile 
model P]. While this model shows many interesting properties, detailed 
investigations § showed that its power spectrum has behaviour in 
any finite dimension. A mean-field calculation of the model did however 
show the expected l//-behaviour exactly p. 

Following the work of BTW a great variety of models (deterministic and 
stochastic, conservative and dissipative, . . . ) have been introduced which 
show the phenomenon of self-organised criticality (SOC). A common feature 
of these models is the presence of a separation of time scales; the system is 
driven at a very slow rate until one of his elements reaches a threshold. This 
triggers a burst of activity (avalanche) which occurs on a very short timescale. 
When the avalanche is over, the system evolves again according to the slow 
drive untill a next avalanche is triggered, .... The activity of the system in 
this way consists of a series of independent avalanches. A generic signature 
of SOC is the presence of a powerlaw in the size (or duration) distribution of 
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the avalanches. If one increases the external driving rate of the system this 
powerlaw disappears. It was however shown by Hwa and Kardar 0, that if 
one increases the rate at which sand is dropped in the sandpile model, and 
one thus obtains the possibility of interacting avalanches, there appears a 
region in the power spectrum where the behaviour is 1//. 

The BTW-sandpile model is a stochastic and conservative model. Olami, 
Feder and Christensen (OFC) [§] introduced a deterministic and dissipative 
model, related to spring-block models of earthquakes, which shows signa- 
tures of SOC, such as the occurence of powerlaw distributions for the sizes 
of the avalanches, with an exponent which depends on the degree of non- 
conservation in the model. In a subsequent study |^ it was shown that this 
model shows 1 / /-noise with an exponent (3 which also depends on the degree 
of non-conservation in the model. In a sense then, the OFC-model fulfills 
more than the sandpile model the original requirements of the concept of 
SOC. 

In the present paper we study the question of l//-noise in the Bak-Sneppen 
model (BS) |]T0|. This model was introduced to describe the coevolution of 
species in the Earth's ecology. Indeed the model shows many qualitative 
similarities with data from the real world, but fails on a quantitative level 



see e.g. \TL\). In this paper we are only interested in the BS-model as an 



interesting physical model and do not discuss its possible biological relevance. 



In the BS-model one has a system of interacting species, each of which 
is represented by a real variable Xi G [0, 1] {i : 1, . . . , A^) which is a mea- 

3 



sure of the fitness of the species. Initially, all Xi are given a random value, 
taken from a uniform distribution on [0, 1]. The dynamics of the model is 
defined as follows. First one looks for the site j where the fitness takes its 
lowest value. One then assigns a new random variable (taken again from the 
uniform distribution) Xj to species j. At the same time, the fitness of K 
other species is changed randomly. Several versions of the BS-model can be 
defined, depending on the way in which these other species are chosen. In 
the lattice version of the model, the species are arranged on a lattice and the 
K species are taken as nearest neighbours. A random neighbour version, in 
which the K neighbours are chosen at random at each timestep, was intro- 
duced in [|12| . This version of the model has the advantage that several of its 



properties can be calculated exactly In this paper we will study both 
this random neighbour version (with K = 1) and a one dimensional version 
of the model in which we only modify the fitness of the neighbour to the 
right of the species with lowest fitness. 

Analytical calculations and extensive simulations have shown that the BS- 
model evolves to a state in which the probability distribution p{x) that a 
species has a fitness x becomes a step function, which is zero for x less then 
some threshold value Xc < 1, and which is 1/(1 — Xc) for x > Xc- In the 
random neighbour model it is known that Xc = 1/{K + 1) exactly. The exact 
value of Xc is not known for any lattice version of the model, but precise 
numerical estimates exist, especially in d = 1, for the case in which both 
neighbours are updated [|n|, |T^]. For the case of the one-dimensional model 
in which one neighbour is updated, we know of no estimate for Xc in the 
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literature. From our numerical results, we estimate Xc = -710 ± 0.005 for 
this case (details of our numerical work are described below). 

Once the system has reached the equilibrium state, its dynamics is charac- 
terised by periods (identified with avalanches) in which at least one of the 
species has a fitness less then Xc , separated by periods in which all species 
have a fitness above threshold. The avalanches can be characterised either 
by their duration or by the their total activity. Let us denote by n{t) the 
number of species which are below threshold as a function of (discrete) time 
t. The total activity s of an avalanche lasting from t = t- io t = tj^ (so its 
total duration is T = t+ — t_ + 1) is then given by 

s = E <t) (1) 

t=t_ 

The distributions P{T) of avalanche durations and P{s) of avalanche sizes 
follow a powerlaw 

P{T) ~ T-^ , P{s) ~ s-y (2) 

For the random neighbour model, it is known exactly that r = 3/2 ||13| while 
for the one dimensional model (2 neighbour updating) the most accurate 
numerical estimate is r = 1.073 ± .003 0]. Our simulations of the one- 
dimensional one neighbour model lead to the estimate r = 1.08 ± .01 giving 
strong evidence that, as could be expected, both one-dimensional models are 
in the same universality class. We don't know of any existing estimates of 
the exponent y for the BS-model. We will return to this exponent at the end 
of the paper. 
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It is of importance to remark that in the BS-model as described so far there 
is no exphcit time separation between a fast time scale for avalanches and 
a slow time scale for inter-avalanche periods. Such a separation is however 
implicitly present in the definition of the model since one assumes that one 
time step in the model is related to a step in 'geological' time tg = exp Xmin/T 
(where Xmin is the lowest value of x at a given time and T is a measure of 
mutation rate, see e.g. |1T|)- When 1/T ^ 1, avalanches occur on time 
scales which are short compared to the timescale of the external drive which 
is set by the mutation rate. 



In order to study spectral properties of the BS-model it is necessary to in- 
troduce a dynamical correlation function 0^(1). In this model a natural 
definition of a correlation function is 

GN{t) = < n{to)n{to + t) >t, - < n{U) >\ (3) 

where the average is taken over time to in the equilibrium state. According 
to the dynamical scaling hypothesis |]T6| one expects the Fouriertransform 
G]^{uj) of a correlation function such as (^) to scale as 

Gm{oj) = uj-^^HiuN') (4) 

where if is a scaling function and z the dynamical exponent. Or equivalently, 
in real space 

Gjv(t) = N'^'^-^^Hit/N') (5) 
We have calculated G'Ar(t) analytically for the random neighbour version 
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{K = 1) of the BS- model and numerically for the one- dimensional one neigh- 
bour version of the model. In both cases we find the presence of l//-noise. 
We now turn to the details of these calculations, and we start with the ana- 
lytical results. 



In a master equation approach to the random neighbour model was 
introduced. Let be the probability that at time t ,n species have a 

fitness which is below a certain value A. In the end we will be most interested 
in the case when X = Xc but for the moment we look at the more general 
case. It is then rather easy to write down a master equation for Pn{t) 

N 

Pn{t + 1) = Mnm Pmit) (6) 

m=0 



where the matrixelements Mnm can be written down in terms of A and [|13 
For t — s> oo, Pn{t) evolves to an equilibrium distribution P*. The correlation 
function G'Ar(t) can also be written down in terms of the matrix M. One has; 

n 2 

(7) 



N N 

G^it) = lim ^ ^mfcP„(to)[M*P(to)]fe 



N 

lim V mPmito) 

to^oo — ' 
m=U 



This expression in fact allows a (numerically) exact calculation of Gi\i{t) 
in finite systems by simple iteration of the master equation (^). We have 
performed such calculations for X = Xc for systems with up to 4000 and 
times t up to 2N (results are discussed below). 

More interesting is the scaling limit in which N ^ oo and A — ^ Xc- In that 
limit it is possible to get a closed expression for the dynamic correlation 
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function. It is therefore convenient to rewrite as 

■ N 



N N 
m=0 k=0 



1 2 



.m=0 



(8) 



where Qmk{t) is the probabihty that in t-timesteps the number of species 



with fitness below A changes from m to k. The authors of assume that 
in the scahng hmit P* becomes a scahng function / of the variable n/\fN 



f 



n 



N , 



(9) 



Inserting (|^) into and taking t ^ oo, — >• oo and A — *• Xc then gives a 
differential equation from which / can be calculated (see eqn. (21) of ||13|| ). 
Using this result we immediately get the second term on the rhs of (Bf) 



N 



* l2 



771=0 



N 
2^ 



(10) 



What remains is a calculation of Qmk{t) in the scaling limit. We therefore 
assume that this probability scales as 



Qmk (^) 



k t 
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(11) 



If we insert this assumption in (^ and take the scaling limit, we obtain a 
differential equation for g (with x = k/\/N,y = m/^/N and r = t/N); 



dg dg 1 d'^g 

dr ^ ^ dx 4 dx"^ 



(12) 



which has to solved with the initial condition 



g{x,y,T = 0) = 6{x - y) 
and reflecting boundary conditions in x = 0. 
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(13) 



The solution is 



9{x, y, t) = h{x, y, r) + h{x, -y, r) 



(14) 



where 



h{x,y,T) 



1 



TT ^1 — exp (— 2r) 
2 



1/2 



exp 



1 - exp {-2t] 



exp 2y'^. 

(y^ + - 2x?/exp (-r))| (15) 



This result has to be used, together with (|TT]), in the first term on the rhs 



of (P). Taking the scaling limit and using the expression of from [0 we 
can rewrite this term as 

.2^2 



TV- 



dx / dy x.y exp {-2y )g{x, y, r) 
Jo 



Inserting our result for g{x, y, r) and performing the integration then finally 
gives; 



GN{t) = iV{^ (1- exp -2r)=^/2[F(l, 2, 3/2, r_(r)) + F(l, 2,3/2, r+(r)) 
-F(l, 2, 5/2, r4r))/3 - F(l, 2, 5/2, r+(r))/3] - ^} (16) 



where F{a, b, c, z) is the hypergeometric function and where 

r±{r) = ^ (1 ±exp-r) 



We thus see that the correlation function has indeed the scaling form (^ with 
z = 1 and a = 2. In figure 1 we show our result (|16D for GN(t)/N versus 
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r, together with the numerical results obtained from direct computation of 
(H) in finite systems. The agreement is perfect thus lending support to the 
scaling assumptions we made. 



To obtain the power spectrum we only have to Fouriertransform (0). Un- 
fortunately, we were not able to obtain an analytical expression for this 
transform. The result of a numerical transform using ©MATHEMATICA 
is shown in figure 2. We show Gn{uj)/N'^ versus uN , which are the natural 
scaling variables according to (^). The straight line shown has a slope — 1. 
These results then show that over many order of magnitude 

Gn{uj) (17) 

SO that indeed there is 1/ /-noise in the model. 

It is interesting to remark here that the random neighbour versions of both 
the BTW-sandpile model and the BS-model [|12| can be related to the 



critical branching process [18|. Within this approximation both models are 
thus in the same universality class. Since it is known that in a mean-field 
theory the sandpile model shows l//-noise it is not so surprising to find 
the same results for the BS-model. 

We now turn to a discussion of the one-dimensional one neighbour version 
of the BS-model. Due to long range correlations which are present between 
subsequent species that have lowest fitness |10| a master equation approach 



is no longer possible. So far, the only approach known for these lattice 
versions of the BS-model is numerical. We have therefore performed extensive 
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numerical calculations of the model on one dimensional lattices with N up 
to 8192 and for time t up to 2^^. Using these data the values of Xc and r for 
the one- dimensional one neighbour model mentionned above were obtained. 
Figure 3 shows numerical results for the correlation function 0^(1) for various 
system sizes. Surprisingly, for large system sizes the correlation function 
seems to become independent of N implying that z becomes 0. We don't 
fully understand this result, but it may be connected with similar behaviour 
found for an other exponent (77) in |[T5|| . 



Figure 4 shows the power spectrum of the correlation function for the system 
with = 8192. As can be seen the behaviour is of the form uj~^ over many 
orders of magnitude. We estimate /? = .97 ± .05. Thus contrary to the 
sandpile model, the BS-model has 1/f behaviour also in a lattice version of 
the model. The exponent f3 is furthermore remarkably close to its mean-field 
value. 

We conclude by interpreting our results in the light of a general theory for 
1 / /-noise for systems with self-organised criticality, put forward in 0,1^. We 
therefore have to introduce first one more exponent, denoted fi which relates 
the (average) duration < T >s of an avalanche to its size s as: 

< T >, ~ sf" (18) 

In it is shown on quite general grounds that when 2fi + r > 3 a model 
shows 1 / /-noise with 

/? = ^ (19) 
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The value of /i for the one- dimensional BS model can be obtained from 
11^ , p!5[] as follows. Let < n(ta) > denote the average number of species below 
threshold a time after the start of an avalanche. In [|14|, |T^ it is shown 
(numerically) that < n{ta) > grows slower then any power (and even becomes 
a constant according to fl^). Therefore, for long living avalanches, it follows 
from (|I|) that < T >s ~ s, or /i = 1 for the one-dimensional model. Turning 
to the random neighbour BS-model, we don't know how to calculate fi exactly 
from the master quation approach to this model. From the precise form of 



the transition probability matrix M as derived in |jT3|, it is however clear that 
for A = Xc, n{t) performs a random walk, from which it can be concluded that 
for long times after the start of an avalanche < n(ta) > ~ t^^^, so that using 
(|l]) we obtain /i = 2/3. We have indeed verfied this /x- value in simulations 
of the random neighbour BS model. (This is an appropriate place to remark 
that these estimates of /i allow us to determine the exponent y (see (0) from 
the obvious relation y = [(r — + 1 so that for the one dimensional 

BS-model we obtain y = t while for the random neighbour model y = 7/4.) 

Using these values we find for the random neighbour version of the BS-model 
2/i + r = 17/6 < 3 so that the results of p predict /? = 2, or absence of 
l//-noise. For the one dimensional model the results of this reference do 
predict l//-noise but with an exponent (3 ~ 1.92. Both predictions are in 
contradiction with the numerical results which we presented here. A possible, 
heuristic, explanation for this is the following. In the arguments of time 
correlations within one avalanche are considered. This is the natural thing to 
do, since in a generic SOC-model the avalanches occur instantly on the long 



12 



time scale, and therefore correlations including the inter-avalanche period 
would be trivially zero. As argued above this explicit time separation is not 
present in the BS-model. Our calculation of the time correlation function 
therefore include periods in which their is no activity. These lead to a more 
rapid decorrelation as compared to the case where one only studies intra- 
avalanche correlations. This more rapid decorrelation in turn leads to a 
decrease in the power of the low frequency components of the power spectrum 
and therefore to a /^-exponent which is lower then that following from simple 
application of the results in 0, . In a sense then, the BS-model is somewhat 
akin to the sandpile model with a finite driving rate as studied in 0, but 
without the possibility of interacting avalanches. This mechanism seems to 
lead to a model which shows both powerlaws in space and time correlations 
without destroying the powerlaws in the size (and/or time) distribution of 
avalanches. It is an interesting project to investigate the generality of this 
scheme. 
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Figure Captions 



Fig.l The exact correlation function of the random neighbour Bak-Sneppen 
model. The figure shows the exact result (16) together with appropriately 
scaled finite system results obtained using (8). The results are for N — 
250(0), N = 500(x), N = lOOO(n), N = 2000(+) and N = 4000(A). 

Fig. 2 Numerical Fouriertransform (open circles) of the exact correlation 
function (12). The straight line represents a best fit through the linear part 
of the data and has a slope of -1. 

Fig. 3 Numerical results for the correlation function of the one-dimensional 
one neighbour Bak-Sneppen model. The different curves represent results for 
(bottom to top) N = 128, 256, 1024 and 4096 respectively. The upper two 
curves almost completely coincide. 

Fig. 4 Fouriertransform (crosses) of the numerically calculated correlation 
function of the one-dimensional one neighbour Bak-Sneppen model. The 
results are for a system oi N — 8192 species. The straight line represents a 
best fit through the hnear part of the data and has a slope of -.972. 
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